library(here)
#analisis
library(ggsignif)
library(ggrepel)
library(ggpubr)
library(inlmisc)
library(nortest) #para testear distribucion
library(skimr) #provides a frictionless approach to summary statistics w
library(easystats) # multiples unciones analiticas
library(lme4)
library(skimr)
library(readxl)
library(fitdistrplus)
# vizualizacion
library(ggridges)
library(sf)
library(GGally)
library(tidyverse, quietly = TRUE)
library(knitr, quietly = TRUE)
library(kableExtra)
library(raster)
library(egg)
library(car) #Variance inflation Factor
library(ggthemes)
library(sjPlot)
library(GGally)Informe autocontenido Análisis Exploratorio de Datos Monitoreo de la pesquería Bacalao de profundidad APA
Analisis exploratorio 1986-2021
ANTECEDENTES
El presente reporte tiene un codigo autocontenido con los Analísis Exploratorios de Datos (AED) de la pesquería de bacalao de profundidad extraído por la flota pesquera artesanal (espinel). Los analisis aca descritos estan compuestos de dos bases. La primera es de los registros oficiales de captura recogidos por el SERNAPESCA. En segundo lugar, se presentan el AED de los datos del monitoreo d la pesqería llevado a cabo por el IFOP a través del Programa de Seguimiento de Pesquerías Demersales del Departamento de Evaluación de Pesquerías.
El principal objetivo es identificar vacios y fortalezas de los datos en terminos espacio-temporales, y con ello tomar decisiones para el paso posterior de modelación de la dinámica poblacional del recurso para realizar asesoría en terminos de estatus y recomendación de una Captura Biologicamente Aceptable (CBA).
ANÁLISIS EXPLORATORIO DE DATOS (AED)
Cargo librerías necesarias para el análisis exploratorio de los datos de las distintas bases de bitácora, tallas y biológico.
# A function for dotplots
multi_dotplot <- function(filename, Xvar, Yvar){
filename %>%
ggplot(aes(x = {{Xvar}})) +
geom_point(aes(y = {{Yvar}}),
alpha=0.4) +
theme_bw() +
coord_flip() +
labs(x = "Order of Data")}Identifico los directorio de trabajo y leo las tres bases de datos; a saber:
- Datos FIP (Rendimiento) 96-32
- Bitácora
- Biologico
- Tallas
- Desembarques (1985-2022)
bit <- read_excel("Artesanal_Historica/datos Bacalao historico 1997-2022CTP/BITACORAS ESPINEL bacalao_corr.xlsx", sheet = "BITACORAS_ESPINEL_bacalao", guess_max = 100000)
long <- read_excel("Artesanal_Historica/datos Bacalao historico 1997-2022CTP/LONGITUD_ESPINEL_bacalao.xlsx")
bio <- read_excel("Artesanal_Historica/datos Bacalao historico 1997-2022CTP/BIOLOGICO ESPINEL bacalao.xlsx",
sheet = "BIOLOGICO_ESPINEL_bacalao")
landing <- read_excel("Artesanal_Historica/DESEMBARQUE HISTÓRICO.xlsx",
skip = 2)
bit8696<- read_csv2("bit_art_86_20.csv")Identificamos los registros asociados recurso objetivo
- Codigo del Recurso: 37
- Nombre común: Bacalao de profundidad
- Nombre científico: Dissostichus eleginoides
Idetifico los registros por recurso y por base para luego filtrar.
dim(bit)
[1] 14208 77
table(bit$COD_ESPECIE)
2 4 5 6 14 15 16 17 18 20 22 23 24
7 1 14 25 147 6 89 16 698 1 25 3 7
25 27 33 35 37 42 43 49 53 56 59 74 75
3 11 2 2 11074 3 9 21 1 5 35 1 28
81 88 99 103 104 105 106 108 119 127 129 144 194
19 1 159 21 24 2 7 44 1 1 1 1 91
199 212 226 239 243 323 324 325 334 844 849 920 984
178 18 1 22 1 2 50 6 830 3 138 1 1
989 999 1309
36 249 58
dim(long)
[1] 15831 24
table(long$COD_ESPECIE)
37
15831
dim(bio)
[1] 127409 36
table(bio$COD_ESPECIE)
37
127409
# las bases de biologicos y longitud solo teienen registros de bacalao.Comenzar a trabajar bases por separado
DESEMBARQUES
Tabla con los desembarques oficiales (Sernapesca, 2022)
kbl(landing, booktabs = T,format = "html",
caption = "Desembarque Bacalao Artesanal por Región") %>%
kable_styling(latex_options = c("striped",
"condensed","scale_down"),
full_width = FALSE) | ...1 | AyP | TPCA | ANTOF | ATCMA | COQ | VALPO | LGBO | MAULE | ÑUBLE | BBIO | ARAUC | RIOS | LAGOS | AYSEN | MAG | ...17 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1985 | NA | 8 | 102 | 146 | 162 | 1674 | NA | 349 | NA | 1527 | NA | NA | 63 | 28 | NA | 4059 |
| 1986 | NA | 842 | 1531 | 213 | 329 | 1242 | NA | 88 | NA | 1546 | NA | NA | 311 | 6 | NA | 6108 |
| 1987 | NA | 182 | 335 | 171 | 51 | 671 | NA | NA | NA | 1264 | NA | NA | 689 | 21 | NA | 3384 |
| 1988 | NA | 131 | 174 | 251 | 42 | 670 | NA | 46 | NA | 1597 | NA | NA | 877 | 8 | NA | 3796 |
| 1989 | NA | 217 | 175 | 177 | 109 | 756 | NA | 114 | NA | 1859 | 10 | NA | 1462 | 8 | NA | 4887 |
| 1990 | NA | 335 | 161 | 126 | 100 | 444 | NA | 5 | NA | 2620 | NA | NA | 1789 | 36 | NA | 5616 |
| 1991 | NA | 280 | 186 | 142 | 223 | 519 | NA | NA | NA | 1901 | NA | NA | 680 | NA | NA | 3931 |
| 1992 | NA | 30 | 45 | 77 | 53 | 477 | NA | 131 | NA | 2108 | 1 | NA | 741 | 1 | NA | 3664 |
| 1993 | NA | 37 | 4 | 66 | 81 | 1676 | NA | 295 | NA | 1111 | NA | NA | 807 | NA | 45 | 4122 |
| 1994 | NA | 24 | 139 | 213 | 68 | 671 | NA | 526 | NA | 1920 | NA | NA | 1374 | 41 | 411 | 5387 |
| 1995 | NA | 252 | 259 | 286 | 65 | 206 | NA | 190 | NA | 1347 | NA | NA | 1680 | 31 | 266 | 4582 |
| 1996 | NA | 476 | 220 | 121 | 177 | 272 | NA | 226 | NA | 1390 | 2 | NA | 1989 | 22 | 83 | 4978 |
| 1997 | NA | 36 | 82 | 47 | 57 | 98 | NA | 116 | NA | 1077 | 1 | NA | 1803 | 22 | 83 | 3422 |
| 1998 | NA | 9 | 67 | 169 | 61 | 155 | NA | 399 | NA | 1284 | NA | NA | 2038 | NA | 11 | 4193 |
| 1999 | NA | 42 | 72 | 168 | 87 | 484 | NA | 406 | NA | 1369 | NA | NA | 2910 | NA | 270 | 5808 |
| 2000 | NA | 285 | 98 | 153 | 105 | 200 | NA | 292 | NA | 800 | NA | NA | 3515 | NA | 345 | 5793 |
| 2001 | NA | 198 | 85 | 127 | 113 | 191 | NA | 136 | NA | 1049 | NA | NA | 2045 | NA | NA | 3944 |
| 2002 | NA | 154 | 44 | 97 | 57 | 188 | NA | 186 | NA | 732 | NA | NA | 3097 | NA | 10 | 4565 |
| 2003 | NA | 58 | 25 | 51 | 27 | 338 | NA | 123 | NA | 752 | NA | NA | 3189 | NA | 179 | 4742 |
| 2004 | NA | 76 | 32 | 63 | 23 | 308 | NA | 132 | NA | 412 | NA | NA | 2303 | NA | 70 | 3419 |
| 2005 | NA | 111 | 14 | 47 | 26 | 84 | NA | 175 | NA | 475 | NA | NA | 2021 | NA | 325 | 3278 |
| 2006 | NA | 87 | 3 | 50 | 50 | 92 | NA | 159 | NA | 243 | NA | NA | 1367 | NA | 40 | 2091 |
| 2007 | 19 | 43 | 7 | 22 | 25 | 43 | NA | 84 | NA | 189 | NA | 501 | 1157 | NA | NA | 2090 |
| 2008 | 4 | 48 | 21 | 9 | 3 | 32 | NA | 70 | NA | 219 | NA | 349 | 803 | NA | NA | 1558 |
| 2009 | 13 | 85 | 24 | 20 | 16 | 79 | NA | 112 | NA | 239 | NA | 522 | 571 | NA | NA | 1681 |
| 2010 | 13 | 43 | 13 | 58 | 11 | 37 | NA | 63 | NA | 109 | NA | 465 | 655 | NA | NA | 1467 |
| 2011 | 19 | 82 | 20 | 46 | 6 | 49 | NA | 90 | NA | 276 | NA | 613 | 987 | NA | NA | 2188 |
| 2012 | 25 | 68 | 11 | 82 | 14 | 26 | NA | 105 | NA | 317 | NA | 290 | 1126 | NA | NA | 2064 |
| 2013 | 26 | 58 | 12 | 68 | 23 | 52 | NA | 214 | NA | 186 | NA | 277 | 622 | NA | 20 | 1558 |
| 2014 | 13 | 48 | - | 59 | 34 | 64 | NA | 171 | NA | 129 | NA | 252 | 453 | NA | 57 | 1280 |
| 2015 | 6 | 52 | 13 | 53 | 38 | 93 | NA | 119 | NA | 261 | NA | 298 | 638 | NA | 38 | 1609 |
| 2016 | 17 | 51 | 20 | 59 | 16 | 70 | NA | 71 | NA | 335 | NA | 327 | 755 | NA | NA | 1721 |
| 2017 | 20 | 92 | 33 | 89 | 33 | 117 | NA | 116 | NA | 316 | NA | 312 | 790 | NA | 28 | 1946 |
| 2018 | 29 | 78 | 5 | 74 | 49 | 78 | NA | 113 | NA | 251 | NA | 218 | 620 | NA | 25 | 1540 |
| 2019 | 26 | 55 | 7 | 63 | 35 | 42 | NA | 102 | NA | 248 | NA | 203 | 696 | NA | 94 | 1571 |
| 2020 | 5 | 18 | 7 | NA | 1 | 27 | NA | 78 | NA | 145 | NA | 229 | 261 | NA | 32 | 803 |
| 2021 | 2 | 83 | 10 | 138 | 17 | 52 | NA | 206 | NA | 304 | NA | 500 | 517 | NA | 21 | 1850 |
| 2022 | 26 | 100 | - | 148 | 58 | 75 | - | 143 | - | 318 | - | 332 | 754 | - | 29 | 1983 |
Ploteo los desembarques por region y por año
landing <- as.data.frame(lapply(landing, as.double))
str(landing)
'data.frame': 38 obs. of 17 variables:
$ ...1 : num 1985 1986 1987 1988 1989 ...
$ AyP : num NA NA NA NA NA NA NA NA NA NA ...
$ TPCA : num 8 842 182 131 217 335 280 30 37 24 ...
$ ANTOF: num 102 1531 335 174 175 ...
$ ATCMA: num 146 213 171 251 177 126 142 77 66 213 ...
$ COQ : num 162 329 51 42 109 100 223 53 81 68 ...
$ VALPO: num 1674 1242 671 670 756 ...
$ LGBO : num NA NA NA NA NA NA NA NA NA NA ...
$ MAULE: num 349 88 NA 46 114 5 NA 131 295 526 ...
$ ÑUBLE: num NA NA NA NA NA NA NA NA NA NA ...
$ BBIO : num 1527 1546 1264 1597 1859 ...
$ ARAUC: num NA NA NA NA 10 NA NA 1 NA NA ...
$ RIOS : num NA NA NA NA NA NA NA NA NA NA ...
$ LAGOS: num 63 311 689 877 1462 ...
$ AYSEN: num 28 6 21 8 8 36 NA 1 NA 41 ...
$ MAG : num NA NA NA NA NA NA NA NA 45 411 ...
$ ...17: num 4059 6108 3384 3796 4887 ...
summary(landing)
...1 AyP TPCA ANTOF
Min. :1985 Min. : 2.00 Min. : 8.00 Min. : 3.00
1st Qu.:1994 1st Qu.:11.25 1st Qu.: 44.25 1st Qu.: 12.75
Median :2004 Median :18.00 Median : 77.00 Median : 32.50
Mean :2004 Mean :16.44 Mean :128.26 Mean : 112.67
3rd Qu.:2013 3rd Qu.:25.25 3rd Qu.:148.25 3rd Qu.: 111.25
Max. :2022 Max. :29.00 Max. :842.00 Max. :1531.00
NA's :22 NA's :2
ATCMA COQ VALPO LGBO
Min. : 9.0 Min. : 1.00 Min. : 26.0 Min. : NA
1st Qu.: 58.0 1st Qu.: 23.50 1st Qu.: 65.5 1st Qu.: NA
Median : 82.0 Median : 49.50 Median : 136.0 Median : NA
Mean :106.7 Mean : 64.34 Mean : 325.1 Mean :NaN
3rd Qu.:148.0 3rd Qu.: 77.75 3rd Qu.: 468.8 3rd Qu.: NA
Max. :286.0 Max. :329.00 Max. :1676.0 Max. : NA
NA's :1 NA's :38
MAULE ÑUBLE BBIO ARAUC RIOS
Min. : 5.0 Min. : NA Min. : 109.0 Min. : 1.0 Min. :203.0
1st Qu.: 99.0 1st Qu.: NA 1st Qu.: 253.5 1st Qu.: 1.0 1st Qu.:270.8
Median :127.0 Median : NA Median : 603.5 Median : 1.5 Median :319.5
Mean :165.3 Mean :NaN Mean : 848.0 Mean : 3.5 Mean :355.5
3rd Qu.:194.0 3rd Qu.: NA 3rd Qu.:1363.5 3rd Qu.: 4.0 3rd Qu.:473.8
Max. :526.0 Max. : NA Max. :2620.0 Max. :10.0 Max. :613.0
NA's :2 NA's :38 NA's :34 NA's :22
LAGOS AYSEN MAG ...17
Min. : 63.0 Min. : 1.00 Min. : 10.00 Min. : 803
1st Qu.: 661.2 1st Qu.: 8.00 1st Qu.: 28.25 1st Qu.:1753
Median : 842.0 Median :22.00 Median : 51.00 Median :3402
Mean :1267.2 Mean :20.36 Mean :112.82 Mean :3228
3rd Qu.:1799.5 3rd Qu.:29.50 3rd Qu.:157.75 3rd Qu.:4472
Max. :3515.0 Max. :41.00 Max. :411.00 Max. :6108
NA's :27 NA's :16
landing2 <- landing %>%
pivot_longer(cols=c("AyP" , "TPCA" ,"ANTOF", "ATCMA", "COQ",
"VALPO" ,"LGBO" , "MAULE", "ÑUBLE", "BBIO" , "ARAUC",
"RIOS" , "LAGOS", "AYSEN", "MAG"),
names_to = "REGION",
values_to = "CAPTURA") %>%
rename(AÑO="...1") %>%
dplyr::select(-2)Genero un gráfico de barras
landing2$REGION <- factor(landing2$REGION ,
levels = c("AyP" , "TPCA" ,"ANTOF", "ATCMA", "COQ",
"VALPO" ,"LGBO" , "MAULE", "ÑUBLE", "BBIO" , "ARAUC",
"RIOS" , "LAGOS", "AYSEN", "MAG"))
desem <- ggplot(landing2 %>%
filter(REGION != "MAG") %>%
drop_na(REGION),aes(AÑO, CAPTURA, fill=REGION)) +
geom_bar(stat="identity")+
scale_fill_viridis_d(option="G")+
theme_few()+
scale_x_continuous(breaks = seq(from = 1985, to = 2022, by = 5))+
scale_y_continuous(breaks = seq(from = 0, to = 4000, by = 1000))+
theme(axis.text.x = element_text(angle = 90, hjust = 2),
panel.grid = element_blank(),
legend.position = "none")+
facet_wrap(~REGION, ncol=7)+
labs(y="Desembarques Oficiales (t)",
x="")
desemLos principales desembarques estan asociados a las regiones de Antifagasta, Valparaiso y BioBIo pero solo durante los primeros años (1985- mediados del 2000). Luego de esto, todas las regiones vieron disminuidos los registros de extracción del recurso.
BITÁCORA
Esta base tiene como principal objetivo obtener un indicador de esfuerzo para calcular un indice de abundancia relativo como la CPUE.
Primero identifico la estructura de la base y filtro el rrecurso bacalao
# filtro bitacoras dejando solo bacalao
bitb <- bit %>%
filter(COD_ESPECIE==37)%>%
mutate(dfp = as.double(dfp),
PESO_corr = as.double(PESO_corr),
prof_med=as.double(prof_med))
dim(bitb)
[1] 11074 77se eliminan 1208-11074 = 3134 registros
glimpse(bitb)ahora las estadísticas descriptivas de las variables de intéres. En este caso dfp, PESO_corr, prof_med
Miramos los oultiers de los datos
# OUTLIERS
#Order data
bitb <- bitb %>%
mutate(order = seq(1:nrow(bitb)))
#Select continuous variables to plot
p1 <- multi_dotplot(bitb, order, dfp)
p2 <- multi_dotplot(bitb, order, PESO_corr)
p3 <- multi_dotplot(bitb, order, prof_med)
p4 <- multi_dotplot(bitb, order, N_TRIPULANTES)
p5 <- multi_dotplot(bitb, order, HORAS_REPOSO)
p6 <- multi_dotplot(bitb, order, NUMERO_DE_ANZUELOS)
#Plot as a grid
ggarrange(p1, p2, p3, p4,p5, p6, nrow = 3)Remuevo outliers
Los datos de profundidades andan bien, Filtro los datos de la variable dfp entre 0 y 30 días. y repito la inspección
bitb2 <- bitb %>%
filter(dfp>0,
dfp<80,
PESO_corr<50000,
N_TRIPULANTES<30,
NUMERO_DE_ANZUELOS<40000)Miro nuenvamente los datos filtrados
# OUTLIERS
#Order data
bitb2 <- bitb2 %>%
mutate(order = seq(1:nrow(bitb2)))
#Select continuous variables to plot
p1a <- multi_dotplot(bitb2, order, dfp)
p2a <- multi_dotplot(bitb2, order, PESO_corr)
p3a <- multi_dotplot(bitb2, order, prof_med)
p4a <- multi_dotplot(bitb2, order, N_TRIPULANTES)
p5a <- multi_dotplot(bitb2, order, HORAS_REPOSO)
p6a <- multi_dotplot(bitb2, order, NUMERO_DE_ANZUELOS)
#Plot as a grid
ggarrange(p1a, p2a, p3a, p4a, p5a, p6a, nrow = 3)Identifico las dimensionesd e la nueva base y la comparo con la anterior
dim(bitb2)
[1] 4745 78
dim(bitb)
[1] 11074 78
table(bitb$año)
1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012
15 557 1046 1072 585 931 576 196 131 83 43 31 57 47 73 20
2013 2014 2015 2016 2017 2018 2019 2020 2021 2022
85 497 506 583 760 1175 474 497 564 470Los registros se minimiza en ciertos años lo cual genera problemas en la fiabilidad de la estimación.
ahora calculamos el rendimiento nominal
CPUE Nominal
Ahora estimo la CPUE, la cual es el rendimoento entre el esfuerzo y la captura. La medida de esfuerzo será dfp. También calculo un rendimiento con la variable de esfuerzo HORAS_REPOSO pero solo como antecedente dado que no tengo muchos datos.
bitb2 <- bitb2 %>%
mutate(CPUE = PESO_corr/dfp) %>%
mutate(CPUE2 =PESO_corr/HORAS_REPOSO)
bitb2$prof_med_cat <- cut(bitb2$prof_med,
breaks = c(500, 1000, 1500,
2000, 3000))# Frequency polygon plot for catch
his <- ggplot(bitb2, aes(CPUE)) +
geom_freqpoly(bins = 5) +
labs(x = "dfp", y = "Frequency") +
theme_bw() +
theme(panel.border = element_rect(colour = "black",
fill=NA, size = 1))+
facet_wrap(~año)+
theme_few()
# Patterns in the variance? (any lack of homogeneity)
point <- ggplot(bitb2 %>%
filter(PESO_corr>1000), aes(x = dfp, y = (PESO_corr))) +
geom_point(shape = 16, size = 5, alpha = 0.6) +
stat_smooth(method = "lm")+
theme(panel.background = element_blank()) +
theme(panel.border = element_rect(fill = NA, size = 1)) +
theme(strip.background = element_rect(fill = "white",
color = "white", size = 1)) +
theme(text = element_text(size=13)) +
xlab("dfp") + ylab("Peso Corregido")
#Plot as a grid
ggarrange(his, point, nrow = 1)Tendencia del esfuerzo dfp a través de los años y por región
effor <- ggplot(bitb2 %>%
group_by(año, REGION_PUERTO_RECALADA) %>%
summarise(DFPM=mean(CPUE)),
aes(año, DFPM))+
geom_point(shape = 16, size = 3, alpha = 0.7)+
scale_x_continuous(breaks = seq(from = 1995, to = 2022, by = 3))+
geom_smooth(method = 'lm',
colour = 'blue',
size = 1.5)+
facet_wrap(~REGION_PUERTO_RECALADA, ncol=5)+
theme_few()+
theme(axis.text.x = element_text(angle = 90, hjust = 2),
panel.grid = element_blank())+
labs(y="Effort (dfp)",
x="")+
ylim(5,500)
efforefformean <- bitb %>%
group_by(año) %>%
summarise(dfpmean= mean(dfp))
write.table(efformean, "dfptable.txt")Tendencia del esfuerzo NUMERO_DE_ANZUELOS a través de los años y por región
efforanz <- ggplot(bitb2 %>%
group_by(año, REGION_PUERTO_RECALADA) %>%
summarise(NANZ=mean(NUMERO_DE_ANZUELOS)),
aes(año, NANZ))+
geom_point(shape = 16, size = 3, alpha = 0.7)+
scale_x_continuous(breaks = seq(from = 1995, to = 2022, by = 1))+
geom_smooth(method = 'lm',
colour = 'green',
size = 1.5)+
facet_wrap(~REGION_PUERTO_RECALADA, ncol=5)+
theme_few()+
theme(axis.text.x = element_text(angle = 90, hjust = 2),
panel.grid = element_blank())+
labs(y="Effort (dfp)",
x="")
efforanzCalculo los ceros de dfp
#CALCULATE NUMBER OF ZEROS
# What is the percentage of zeros i the response variable
round(sum(bitb2$dfp == 0) * 100 / nrow(bitb2),0)
[1] 0
#0Interacciones
# Interactions
# Year x season
ggplot(bitb2, aes(x = dfp, y = (PESO_corr))) +
geom_point(shape = 16, size = 3, alpha = 0.7) +
geom_smooth(method = 'lm', colour = 'red', se = FALSE) +
theme_bw() +
xlab("Points") + ylab("Catch") +
facet_grid(año~trim)# No
ggplot(bitb2, aes(x = dfp, y = (PESO_corr))) +
geom_point(shape = 16, size = 3, alpha = 0.7) +
geom_smooth(method = 'lm', colour = 'red', se = FALSE) +
theme_bw() +
xlab("dfp") + ylab("Catch") +
facet_grid(año~N_TRIPULANTES)+
theme_few()# Perhaps
ggplot(bitb2, aes(x = dfp, y = (PESO_corr))) +
geom_point(shape = 16, size = 3, alpha = 0.7) +
geom_smooth(method = 'lm', colour = 'red', se = FALSE) +
theme_bw() +
xlab("dfp") + ylab("Catch") +
facet_grid(trim~N_TRIPULANTES)+
theme_few()# CPUE slope varies between habitats - interactionProfundidad por region
profu <- ggplot(bitb2, aes(año, desc(prof_med),
group=año))+
geom_boxplot(fill=NA, alpha=.5)+
geom_jitter(size=0.4, alpha=0.2,
width = .25)+
facet_wrap(~REGION_PUERTO_RECALADA, ncol=5)+
theme_few()+
scale_x_continuous(breaks = seq(from = 2011, to = 2022, by = 3))+
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
labs(x="Región",
y="Profundidad (mts)")
profuComparo las medias de la profundidad
# Prueba de Kruskal-Wallis
pairwise.t.test(x = bitb2$prof_med,
g = bitb2$REGION_PUERTO_RECALADA,
p.adjust.method = "holm",
pool.sd = TRUE,
paired = FALSE,
alternative = "two.sided")
Pairwise comparisons using t tests with pooled SD
data: bitb2$prof_med and bitb2$REGION_PUERTO_RECALADA
1 3 4 5 7 8 10 12 14
3 1.00000 - - - - - - - -
4 1.00000 1.00000 - - - - - - -
5 3.2e-07 0.09875 1.00000 - - - - - -
7 1.00000 1.00000 1.00000 8.0e-11 - - - - -
8 5.5e-08 0.02686 1.00000 1.00000 5.2e-11 - - - -
10 0.00032 0.43592 1.00000 1.00000 0.00012 1.00000 - - -
12 0.68271 0.27691 0.21334 0.01917 0.44715 0.01135 0.02686 - -
14 1.00000 1.00000 1.00000 0.00156 1.00000 0.00027 0.04753 0.35077 -
15 2.1e-12 1.2e-13 1.7e-07 < 2e-16 < 2e-16 < 2e-16 < 2e-16 1.00000 9.3e-15
P value adjustment method: holmggplot(bitb2,
aes(x = nombre_puerto_recalada ,
y = dfp,
fill = nombre_puerto_recalada )) +
geom_violindot(fill_dots = "black") +
theme_modern() +
scale_fill_material_d()
ggplot(bitb2, aes(x = as.factor(REGION_PUERTO_RECALADA ),
y = bitb2$NUMERO_DE_ANZUELOS,
fill = as.factor(REGION_PUERTO_RECALADA))) +
geom_violin(scale="width") +
theme_modern() +
scale_fill_material_d(palette="ice",
name="REGION")ggplot(bitb2 %>%
drop_na(prof_med_cat),
aes(x = dfp, y = (PESO_corr))) +
geom_point(shape = 16, size = 3, alpha = 0.7) +
geom_smooth(method = 'lm',
colour = 'red',
se = FALSE) +
theme_few() +
xlab("dfp") + ylab("Catch") +
facet_grid(prof_med_cat~año)ploteo los datos de CPUE totales.
cpuenom <- ggplot(bitb2 %>%
group_by(zona, año) %>%
summarise(CPUEM=mean(CPUE)) %>%
filter(CPUEM<600,
año>1996,
zona!=4),
aes(año, CPUEM, color=zona, group=zona))+
geom_point()+
geom_smooth()+
scale_color_viridis_b(option="B")+
scale_x_continuous(breaks = seq(from = 2010, to = 2022, by = 1))+
theme_few()+
labs(x="",
y="CPUE")
cpuenomguardo los datos bitacora CPUE nominal 1996-2022 del objeto cpuenom
write_csv(cpuenom, "CPUE_MON.csv")CPUE Estandarizada
Identificar los proncipales factores para modelar la variable.
Luego de una reunión con Seguimiento (Patricio Galvez) se indicaron algunos aspectos que deben ser considerados en la estandariación.
El poder de pesca (tipo de embarcación) es mayor en la zona centro sur de Chile. Identificar equilibrio de la base
bitb2respecto a este dato.Usar zonas como factor en desmedro de región.
Identifico una estadistica descriptiva general de la base de estandarización
Correlaciones
Análisis de utilidad para identificar a traves de un metodo de corrrelación de pearson, la correlación en tre las variables que serán utilizadas en la estandariacion de la CPUE.
Primero identifico las variables a correlacionar
names(bitb2)
[1] "ID" "embarque"
[3] "COD_BARCO" "FECHA_HORA_RECALADA"
[5] "FECHA_HORA_ZARPE" "COD_PESQUERIA"
[7] "COD_PESQUERIA_ant" "N_TRIPULANTES"
[9] "N_CALADAS" "PUERTO_ZARPE"
[11] "PUERTO_RECALADA" "nombre_puerto_recalada"
[13] "NRO_FORMULARIO" "REGION_PUERTO_RECALADA"
[15] "ESPECIE_OBJETIVO" "ECOSONDA"
[17] "GPS" "POTENCIA_MOTOR_EQ"
[19] "VIRADOR" "GASTOS_VIAJE"
[21] "NUMERO_LANCE_EX" "FECHA_LANCE"
[23] "año" "mes"
[25] "trim" "dfp"
[27] "ID_CUADRICULA" "ID_PROCEDENCIA"
[29] "PESO_TOTAL_CAPTURA" "LATITUD"
[31] "LONGITUD" "lat_ctgr"
[33] "lon_ctgr" "zona"
[35] "ESPECIE_OBJETIVO_LANCE" "CLASE_LANCE"
[37] "COD_ESPECIE" "PESO"
[39] "PESO_corr" "CAPTURA_RETENIDA"
[41] "VOLUMEN_DESCARTE" "LUGAR_DESCARTE"
[43] "TIPO_DESCARTE" "PRECIO_UNITARIO"
[45] "DESTINO_CAPTURA" "NRO_ANZUELOS"
[47] "REGION_PROCEDENCIA" "FECHA_HORA_FIN_CALADO"
[49] "FECHA_HORA_INI_CALADO" "LATITUD_INICIO_CALADO"
[51] "LONGITUD_INICIO_CALADO" "LATITUD_FIN_CALADO"
[53] "LONGITUD_FIN_CALADO" "FECHA_HORA_INI_VIRADO"
[55] "FECHA_HORA_FIN_VIRADO" "LATITUD_INICIAL_VIRADO"
[57] "LONGITUD_INICIAL_VIRADO" "LATITUD_FIN_VIRADO"
[59] "LONGITUD_FIN_VIRADO" "PROFUNDIDAD_MINIMA_LM"
[61] "PROFUNDIDAD_MAXIMA_LM" "PROFUNDIDAD_FONDO_VIRADO_INI"
[63] "PROFUNDIDAD_FONDO_VIRADO_FIN" "prof_med"
[65] "SEPARACION_ANZUELO" "MARCA_ANZUELO"
[67] "NUMERO_DE_ANZUELOS" "PERDIDA_ANZUELOS"
[69] "TIPO_CARNADA" "TAMANIO_ANZUELO"
[71] "INTENSIDAD_VIENTO" "LONGITUD_LINEA_MADRE"
[73] "HORAS_REPOSO" "TIPO_ESPINEL"
[75] "N_PANOS" "ANZ_POR_PANO"
[77] "filtro" "order"
[79] "CPUE" "CPUE2"
[81] "prof_med_cat"
bitb3 <- bitb2 %>%
dplyr::select(8, 18, 23,24, 25, 26, 32, 33, 34, 39, 69, 73, 75, 78, 79, 80, 81)
names(bitb3)
[1] "N_TRIPULANTES" "POTENCIA_MOTOR_EQ" "año"
[4] "mes" "trim" "dfp"
[7] "lat_ctgr" "lon_ctgr" "zona"
[10] "PESO_corr" "TIPO_CARNADA" "HORAS_REPOSO"
[13] "N_PANOS" "order" "CPUE"
[16] "CPUE2" "prof_med_cat"results <- correlation(bitb3)
results %>%
summary(redundant=TRUE) %>%
plot(results, show_data = "points")+
theme_bw()skimr::skim(bitb3)| Name | bitb3 |
| Number of rows | 4745 |
| Number of columns | 17 |
| _______________________ | |
| Column type frequency: | |
| character | 1 |
| factor | 1 |
| numeric | 15 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| TIPO_CARNADA | 122 | 0.97 | 1 | 5 | 0 | 8 | 0 |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| prof_med_cat | 939 | 0.8 | FALSE | 4 | (1e: 1867, (50: 961, (1.: 853, (2e: 125 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| N_TRIPULANTES | 0 | 1.00 | 7.23 | 1.62 | 2.00 | 6.00 | 7.00 | 9.00 | 10.00 | ▁▂▆▇▆ |
| POTENCIA_MOTOR_EQ | 167 | 0.96 | 405.74 | 399.03 | 0.00 | 320.00 | 390.00 | 420.00 | 3903.00 | ▇▁▁▁▁ |
| año | 0 | 1.00 | 2018.05 | 2.47 | 2011.00 | 2016.00 | 2018.00 | 2020.00 | 2022.00 | ▁▅▆▆▇ |
| mes | 0 | 1.00 | 7.29 | 3.79 | 1.00 | 3.00 | 9.00 | 11.00 | 12.00 | ▅▂▁▂▇ |
| trim | 0 | 1.00 | 2.71 | 1.32 | 1.00 | 1.00 | 3.00 | 4.00 | 4.00 | ▆▂▁▂▇ |
| dfp | 0 | 1.00 | 16.73 | 6.62 | 0.67 | 13.02 | 16.79 | 19.79 | 50.67 | ▂▇▂▁▁ |
| lat_ctgr | 51 | 0.99 | -35.49 | 6.89 | -55.17 | -40.99 | -35.25 | -33.17 | -18.33 | ▁▅▇▂▂ |
| lon_ctgr | 51 | 0.99 | -73.24 | 1.62 | -81.89 | -74.69 | -73.08 | -71.97 | -51.94 | ▁▇▁▁▁ |
| zona | 0 | 1.00 | 2.11 | 0.65 | 1.00 | 2.00 | 2.00 | 2.00 | 4.00 | ▂▇▁▃▁ |
| PESO_corr | 0 | 1.00 | 755.69 | 1304.56 | 0.00 | 41.00 | 110.00 | 1079.66 | 14092.00 | ▇▁▁▁▁ |
| HORAS_REPOSO | 2988 | 0.37 | 1629.72 | 832.08 | 12.00 | 1200.00 | 1200.00 | 2400.00 | 7200.00 | ▇▅▁▁▁ |
| N_PANOS | 1833 | 0.61 | 226.33 | 280.51 | 1.00 | 1.00 | 2.00 | 420.00 | 900.00 | ▇▁▃▁▁ |
| order | 0 | 1.00 | 2373.00 | 1369.91 | 1.00 | 1187.00 | 2373.00 | 3559.00 | 4745.00 | ▇▇▇▇▇ |
| CPUE | 0 | 1.00 | 54.95 | 92.65 | 0.00 | 2.44 | 6.49 | 84.01 | 1067.68 | ▇▁▁▁▁ |
| CPUE2 | 2988 | 0.37 | 1.55 | 5.23 | 0.00 | 0.51 | 0.99 | 1.79 | 158.44 | ▇▁▁▁▁ |
Datos 0
Are there missing values? y sirve para identificar los factores que se utilzaran en la estandarización.
dim(bitb3)
[1] 4745 17
colSums(is.na(bitb3))
N_TRIPULANTES POTENCIA_MOTOR_EQ año mes
0 167 0 0
trim dfp lat_ctgr lon_ctgr
0 0 51 51
zona PESO_corr TIPO_CARNADA HORAS_REPOSO
0 0 122 2988
N_PANOS order CPUE CPUE2
1833 0 0 2988
prof_med_cat
939Es probable que HORAS_REPOSOde la red no podamos usarlo como factor ni como variable de esfuerzo por lo exiguo del dato.
Variabilidad
asd1<-ggplot(data = bitb3 %>%
filter(zona!=4),
aes(y=dfp, x=zona,
group=zona)) +
geom_boxplot()+
labs(title = "Variabilidad en DFP",
y = "DFP", x = "ZONA") +
facet_wrap(~ año, scales = "free")+
theme_few()
asd1asd2<-ggplot(data =bitb3,
aes(x=LAT_CAT,
y=CPUE)) +
geom_boxplot()+
labs(title = "Variabilidad en distribución latitudinal",
y = "CPUE", x = "Zona") +
theme_few()+
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
facet_wrap(~ año)
asd2Normality and Homogenity factors
# Frequency polygon plot for catch
ggplot(bitb3 %>%
filter(zona!=4),
aes(CPUE)) +
geom_freqpoly(bins = 8) +
labs(x = "Cpue", y = "Frequency") +
facet_wrap(.~zona)+
theme_few()+
theme(panel.border =
element_rect(colour = "black",
fill=NA, size = 1))Shapiro-Wilk test for deviation from normality
shapiro.test(bitb3$CPUE)
Shapiro-Wilk normality test
data: bitb3$CPUE
W = 0.63039, p-value < 2.2e-16Departure from normality…but not critical
Colinealidad entre factores
# COLLINEARITY
Coll <- c("año", "dfp" ,"lat_ctgr", "zona" ,
"N_PANOS" , "CPUE", "prof_med_cat")
# Obtain summary using the ggpairs command from the GGally library
ggpairs(bitb3[,Coll], ggplot2::aes(alpha = 0.9, colour=zona),
lower = list(combo = "count"))+
scale_fill_manual(values=c("red", "blue", "green", "black")) +
scale_colour_manual(values=c("red", "blue", "green", "black")) +
theme_few()
# Nothing serious#Calculate Variance Inflation Factor (VIF)
vifmodel <- round(vif(lm(PESO_corr ~ zona + dfp + año + mes,
data = bitb3)),2)
vifmodel
zona dfp año mes
1.34 1.29 1.06 1.01barplot(vifmodel, main = "VIF Values", horiz = TRUE, col = "steelblue")La multicolinealidad en el análisis de regresión se produce cuando dos o más variables predictoras están altamente correlacionadas entre sí, de modo que no proporcionan información única o independiente en el modelo de regresión.
Si el grado de correlación es lo suficientemente alto entre las variables, puede causar problemas al ajustar e interpretar el modelo de regresión.
La forma más común de detectar la multicolinealidad es mediante el factor de inflación de varianza (VIF), que mide la correlación y la fuerza de la correlación entre las variables predictoras en un modelo de regresión.
El valor de VIF comienza en 1 y no tiene límite superior. Una regla general para interpretar VIF es la siguiente:
Un valor de 1 indica que no hay correlación entre una variable predictora dada y cualquier otra variable predictora en el modelo. Un valor entre 1 y 5 indica una correlación moderada entre una variable predictora dada y otras variables predictoras en el modelo, pero esto a menudo no es lo suficientemente grave como para requerir atención. Un valor superior a 5 indica una correlación potencialmente grave entre una variable predictora dada y otras variables predictoras en el modelo. En este caso, es probable que las estimaciones de coeficientes y los valores p en el resultado de la regresión no sean fiables.
Interaccion
# Patterns in the variance? (Evidence for lack of homogeneity)
ggplot(bitb3,
aes(x = dfp,
y = (CPUE))) +
geom_jitter(shape = 16,
size = 2.5,
alpha = 0.3,
height = 0.25,
width = 0.5) +
geom_smooth(method = 'lm', colour = 'red', se = FALSE, size = 1.5) +
facet_grid(bitb3$N_TRIPULANTES~bitb3$zona)+
theme(panel.background = element_blank()) +
theme(panel.border = element_rect(fill = NA,
size = 1)) +
theme(strip.background =
element_rect(fill = "white",
color = "white",
size = 1)) +
theme(text = element_text(size=13)) +
theme_few()+
xlab("Effort") +
ylab("CPUE") +
ggtitle("variabless por Numero de Tripulantes y por zona")Year x mes
ggplot(bitb3, aes(x = dfp, y = (PESO_corr))) +
geom_point(shape = 16, size = 3, alpha = 0.7) +
geom_smooth(method = 'lm', colour = 'red', se = FALSE, size = 1.5) +
theme_few() +
xlab("dfp") + ylab("Catch") +
facet_grid(año~trim)# No# Zona
ggplot(bitb3, aes(x = dfp, y = (PESO_corr))) +
geom_point(shape = 16, size = 3, alpha = 0.7) +
geom_smooth(method = 'lm', colour = 'red', se = FALSE, size = 1.5) +
theme_few() +
xlab("dfp") + ylab("Catch") +
facet_grid(año~trim) facet_grid(~zona)
<ggproto object: Class FacetGrid, Facet, gg>
compute_layout: function
draw_back: function
draw_front: function
draw_labels: function
draw_panels: function
finish_data: function
init_scales: function
map_data: function
params: list
setup_data: function
setup_params: function
shrink: TRUE
train_scales: function
vars: function
super: <ggproto object: Class FacetGrid, Facet, gg>
# CPUE slope varies between habitats - interaction# Tipo carnada zona
ggplot(bitb3, aes(x = dfp, y = (PESO_corr))) +
geom_point(shape = 16, size = 3, alpha = 0.7) +
geom_smooth(method = 'lm', colour = 'red', se = FALSE, size = 1.5) +
theme_few()+
xlab("Points") + ylab("Catch") +
facet_grid(TIPO_CARNADA~zona)# PerhapsGenero categorias y Factorizo
bitb3 <- bitb3 %>%
mutate(POTENCIA_CAT = cut(POTENCIA_MOTOR_EQ,
breaks = c(0, 200, 300, 400, 1000),
labels = c("Bajo", "Medio-Bajo", "Medio-Alto", "Alto")),
LAT_CAT = cut(lat_ctgr,
breaks = c(-20, -35, -40, -47),
labels = c("Norte", "Centro", "Sur")),
N_TRIPULANTES=factor(N_TRIPULANTES),
año=factor(año),
zona=factor(zona),
trim=factor(trim)) %>%
drop_na(zona,
LAT_CAT,
prof_med_cat,
POTENCIA_CAT) %>%
filter(CPUE<600) %>%
mutate(logCPUE=log(CPUE))
glimpse(bitb3)
Rows: 3,340
Columns: 20
$ N_TRIPULANTES <fct> 4, 8, 8, 8, 6, 7, 6, 7, 8, 8, 6, 8, 8, 9, 9, 9, 9, 9…
$ POTENCIA_MOTOR_EQ <dbl> 250, 380, 320, 450, 380, 420, 370, 240, 320, 360, 45…
$ año <fct> 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013…
$ mes <dbl> 7, 7, 7, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 11,…
$ trim <fct> 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4…
$ dfp <dbl> 13.833333, 11.741667, 15.086806, 13.729167, 18.31736…
$ lat_ctgr <dbl> -40.66667, -39.96667, -40.50000, -33.91667, -40.6666…
$ lon_ctgr <dbl> -74.50000, -74.33333, -74.50000, -72.58333, -74.6666…
$ zona <fct> 2, 2, 2, 2, 2, 3, 3, 3, 3, 2, 3, 3, 3, 2, 2, 2, 2, 2…
$ PESO_corr <dbl> 1095.0568, 1256.6812, 1422.8584, 975.6000, 2026.6464…
$ TIPO_CARNADA <chr> "2", "3", "3", NA, "1", "2", "2", "2", "3", "3", "3"…
$ HORAS_REPOSO <dbl> 1400, 1800, 1200, 4800, 1200, 1200, 2400, 1200, 1200…
$ N_PANOS <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ order <int> 8, 10, 12, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 2…
$ CPUE <dbl> 79.160730, 107.027493, 94.311443, 71.060393, 110.640…
$ CPUE2 <dbl> 0.7821834, 0.6981562, 1.1857153, 0.2032500, 1.688872…
$ prof_med_cat <fct> "(1e+03,1.5e+03]", "(500,1e+03]", "(1e+03,1.5e+03]",…
$ POTENCIA_CAT <fct> Medio-Bajo, Medio-Alto, Medio-Alto, Alto, Medio-Alto…
$ LAT_CAT <fct> Norte, Centro, Norte, Sur, Norte, Norte, Norte, Nort…
$ logCPUE <dbl> 4.371480, 4.673086, 4.546603, 4.263530, 4.706288, 4.…The data exploration showed:
- One significant outlier in catch
- Deviation from normality in response variable
- Possible departure from homogeneity
- Few zeros in the response variable
- No collinearity
- No imbalance
- Possible season x habitat interaction
Chequeo distribuciones
bcpue <- ggplot(bitb3,
aes(CPUE, group=año))+
coord_flip()+
geom_boxplot()+
theme_few()
hcpue <- ggplot(bitb3 ,
aes(CPUE))+
geom_histogram(bins=15, fill=2)+
theme_few()
ggarrange(bcpue, hcpue, ncol = 2)Grafico para comprobar distribución Gamma
dist <- descdist(bitb3$CPUE, boot = 200)Gamma distribution es la mas adecuada para la variable y que tiene la forma de :
f(x, u, \sigma) = \frac{1}{\sqrt{2hr}}e^{\frac{1(x-u)^2}{\sigma^2}}
Donde;
x \in R
En la transformación de la dsata se comprueba la inconcistencia del uso de datos de la variable CPUE transformada
Modelos GLM
#Predictores (Factores Principales)
M01 = glm(formula = CPUE ~ año,
family = gaussian(link = "identity"),
data = bitb3,na.action=na.exclude)
M02 = glm(formula= CPUE ~ año+trim,
family = gaussian(link = "identity"),
data = bitb3,na.action=na.exclude)
M03 = glm(formula= CPUE ~ año+trim+prof_med_cat,
family = gaussian(link = "identity"),
data = bitb3,na.action=na.exclude)
M04 = glm(formula= CPUE ~ año+trim+prof_med_cat+N_TRIPULANTES,
family = gaussian(link = "identity"),
data = bitb3,na.action=na.exclude)
M05 = glm(formula= CPUE ~ año+trim+prof_med_cat+N_TRIPULANTES+LAT_CAT,
family = gaussian(link = "identity"),
data = bitb3,na.action=na.exclude)
M06 = glm(formula= CPUE ~ año+trim+prof_med_cat+N_TRIPULANTES+LAT_CAT+TIPO_CARNADA,
family = gaussian(link = "identity"),
data = bitb3,na.action=na.exclude)
M07 = glm(formula= CPUE ~ año+trim+prof_med_cat+N_TRIPULANTES+LAT_CAT+zona,
family = gaussian(link = "identity"),
data = bitb3,na.action=na.exclude)
M08 = glm(formula= CPUE ~ año+trim+prof_med_cat+N_TRIPULANTES+LAT_CAT+zona:prof_med_cat,
family = gaussian(link = "identity"),
data = bitb3,na.action=na.exclude)rmodelo01 <- c(AIC(M01),(M01$null.deviance-M01$deviance)/M01$null.deviance)
rmodelo02 <- c(AIC(M02),(M02$null.deviance-M02$deviance)/M02$null.deviance)
rmodelo03 <- c(AIC(M03),(M03$null.deviance-M03$deviance)/M03$null.deviance)
rmodelo04 <- c(AIC(M04),(M04$null.deviance-M04$deviance)/M04$null.deviance)
rmodelo05 <- c(AIC(M05),(M05$null.deviance-M05$deviance)/M05$null.deviance)
rmodelo06 <- c(AIC(M06),(M06$null.deviance-M06$deviance)/M06$null.deviance)
rmodelo07 <- c(AIC(M07),(M07$null.deviance-M07$deviance)/M07$null.deviance)
rmodelo08 <- c(AIC(M08),(M08$null.deviance-M08$deviance)/M08$null.deviance)resultados <- as.data.frame(rbind(rmodelo01,rmodelo02,rmodelo03,rmodelo04,rmodelo05,rmodelo06,
rmodelo07,rmodelo08))
resultados <- resultados %>%
rename("AIC"=V1,
"Deviance"=V2)
kbl(resultados, booktabs = T,format = "html",
caption = "Parámetros de los modelos estimados") %>%
kable_styling(latex_options = c("striped",
"condensed","scale_down"),
full_width = FALSE) predict.3 <- predict(M07,se.fit=T)
predict.mean3 <-exp(predict.3$fit)
predict.var3 <-predict.mean3^2*predict.3$se.fit^2
predict.sd3 <-sqrt(predict.var3)
predict.cv3 <-predict.sd3/predict.mean3
predict.linf3 <-predict.mean3-1.96*predict.sd3
predict.lsup3 <-predict.mean3+1.96*predict.sd3
hist(exp(predM06$fit))Seleccion del modelo final M06
check_model(M07)Table with significance level (***).
tab_model(M01,
M02,
M03,
M04,
M05,
M06,
M07,
M08,
p.style = "stars")| CPUE | CPUE | CPUE | CPUE | CPUE | CPUE | CPUE | CPUE | |||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | CI | Estimates | CI | Estimates | CI | Estimates | CI | Estimates | CI | Estimates | CI | Estimates | CI | Estimates | CI |
| (Intercept) | 45.27 *** | 22.46 – 68.08 | 42.90 *** | 19.14 – 66.66 | 33.63 ** | 10.36 – 56.90 | -3.00 | -60.39 – 54.39 | 31.30 | -25.96 – 88.56 | 121.42 | -44.67 – 287.50 | 30.55 | -27.52 – 88.61 | -4.34 | -67.24 – 58.55 |
| año [2014] | -19.05 | -44.15 – 6.04 | -24.71 | -50.17 – 0.76 | -23.76 | -48.36 – 0.83 | -17.55 | -42.00 – 6.90 | -1.95 | -26.33 – 22.43 | -45.62 * | -87.06 – -4.19 | -3.96 | -28.43 – 20.51 | -8.07 | -32.61 – 16.47 |
| año [2015] | -4.10 | -29.54 – 21.34 | -4.61 | -30.14 – 20.92 | -17.57 | -42.24 – 7.10 | -0.00 | -24.87 – 24.87 | 14.90 | -9.94 – 39.73 | -25.45 | -66.99 – 16.09 | 10.86 | -14.18 – 35.90 | 6.08 | -19.01 – 31.17 |
| año [2016] | 31.01 * | 5.22 – 56.80 | 22.00 | -4.86 – 48.87 | 4.41 | -21.56 – 30.39 | 11.85 | -13.89 – 37.60 | 26.23 * | 0.63 – 51.83 | -18.39 | -60.45 – 23.67 | 21.21 | -4.63 – 47.05 | 18.48 | -7.27 – 44.23 |
| año [2017] | 6.25 | -17.66 – 30.17 | 6.30 | -17.68 – 30.28 | -4.39 | -27.54 – 18.77 | -0.40 | -23.49 – 22.68 | 10.37 | -12.54 – 33.27 | -37.39 | -77.88 – 3.11 | 6.12 | -17.10 – 29.34 | 2.09 | -21.17 – 25.35 |
| año [2018] | -3.95 | -27.68 – 19.78 | -4.09 | -27.86 – 19.69 | -23.24 * | -46.28 – -0.20 | -13.79 | -36.88 – 9.30 | -2.77 | -25.70 – 20.16 | -53.80 ** | -94.29 – -13.32 | -10.32 | -33.51 – 12.87 | -15.09 | -38.32 – 8.13 |
| año [2019] | 21.67 | -2.77 – 46.12 | 21.45 | -2.98 – 45.88 | 9.17 | -14.42 – 32.77 | 17.39 | -6.40 – 41.19 | 25.62 * | 2.05 – 49.19 | -25.36 | -66.17 – 15.46 | 20.48 | -3.08 – 44.05 | 14.89 | -8.80 – 38.58 |
| año [2020] | -21.96 | -46.10 – 2.19 | -20.98 | -45.22 – 3.25 | -29.77 * | -53.14 – -6.40 | -11.82 | -35.41 – 11.77 | 3.14 | -20.37 – 26.65 | -48.70 * | -89.56 – -7.83 | -2.38 | -26.16 – 21.40 | -5.47 | -29.31 – 18.37 |
| año [2021] | 55.29 *** | 31.14 – 79.43 | 49.12 *** | 24.71 – 73.52 | 40.88 *** | 17.34 – 64.42 | 48.56 *** | 24.94 – 72.18 | 58.81 *** | 34.95 – 82.68 | 8.87 | -31.94 – 49.67 | 48.70 *** | 24.62 – 72.78 | 48.36 *** | 24.23 – 72.49 |
| año [2022] | 33.32 ** | 8.43 – 58.21 | 27.27 * | 1.95 – 52.58 | 20.30 | -4.10 – 44.71 | 29.96 * | 5.51 – 54.40 | 40.26 ** | 15.55 – 64.97 | -1.15 | -42.46 – 40.16 | 31.55 * | 6.67 – 56.44 | 29.76 * | 4.87 – 54.66 |
| trim [2] | 15.86 ** | 6.07 – 25.65 | 8.40 | -1.08 – 17.88 | 11.33 * | 1.79 – 20.87 | 11.17 * | 1.75 – 20.60 | 11.09 * | 1.76 – 20.43 | 11.02 * | 1.63 – 20.41 | 10.47 * | 1.12 – 19.82 | ||
| trim [3] | 19.61 ** | 7.44 – 31.78 | 16.79 ** | 5.06 – 28.53 | 21.93 *** | 10.19 – 33.67 | 18.50 ** | 6.89 – 30.12 | 14.22 * | 2.70 – 25.74 | 17.35 ** | 5.77 – 28.93 | 17.10 ** | 5.51 – 28.70 | ||
| trim [4] | 1.31 | -5.71 – 8.34 | -2.33 | -9.12 – 4.46 | -1.17 | -7.96 – 5.62 | -1.56 | -8.26 – 5.14 | -4.48 | -11.12 – 2.17 | -0.57 | -7.27 – 6.12 | 0.45 | -6.26 – 7.17 | ||
| prof med cat [1001-1.5e+03] |
17.80 *** | 11.16 – 24.44 | 22.47 *** | 15.72 – 29.22 | 21.81 *** | 15.14 – 28.48 | 23.82 *** | 17.15 – 30.50 | 20.27 *** | 13.60 – 26.95 | 82.41 *** | 55.05 – 109.78 | ||||
| prof med cat [1501-2e+03] | 63.66 *** | 55.44 – 71.88 | 63.17 *** | 54.86 – 71.49 | 59.62 *** | 51.33 – 67.91 | 60.37 *** | 52.14 – 68.61 | 55.48 *** | 47.11 – 63.85 | 74.71 *** | 47.06 – 102.36 | ||||
| prof med cat [2001-3e+03] | 66.83 *** | 47.76 – 85.90 | 68.40 *** | 49.52 – 87.27 | 65.25 *** | 46.60 – 83.90 | 65.75 *** | 47.33 – 84.17 | 59.33 *** | 40.65 – 78.01 | 64.31 ** | 21.80 – 106.82 | ||||
| N TRIPULANTES [4] | 13.90 | -41.73 – 69.54 | 14.88 | -40.06 – 69.82 | 17.14 | -37.06 – 71.34 | 15.78 | -38.89 – 70.45 | 15.63 | -38.74 – 70.01 | ||||||
| N TRIPULANTES [5] | 60.36 * | 6.97 – 113.75 | 54.09 * | 1.34 – 106.85 | 58.17 * | 6.17 – 110.18 | 59.18 * | 6.63 – 111.72 | 58.84 * | 6.52 – 111.16 | ||||||
| N TRIPULANTES [6] | 26.94 | -25.53 – 79.42 | 21.62 | -30.31 – 73.54 | 29.20 | -22.01 – 80.41 | 38.35 | -13.76 – 90.47 | 40.92 | -11.10 – 92.93 | ||||||
| N TRIPULANTES [7] | 5.36 | -47.15 – 57.86 | -4.03 | -56.03 – 47.96 | 4.67 | -46.62 – 55.96 | 17.25 | -35.16 – 69.66 | 19.40 | -32.88 – 71.68 | ||||||
| N TRIPULANTES [8] | 31.78 | -20.95 – 84.51 | 12.54 | -39.88 – 64.96 | 16.01 | -35.68 – 67.69 | 31.99 | -20.74 – 84.72 | 33.38 | -19.22 – 85.98 | ||||||
| N TRIPULANTES [9] | 37.40 | -15.52 – 90.33 | 3.70 | -49.26 – 56.66 | 14.32 | -37.94 – 66.58 | 21.57 | -31.54 – 74.68 | 23.66 | -29.29 – 76.60 | ||||||
| N TRIPULANTES [10] | 1.39 | -56.22 – 59.00 | -39.17 | -96.77 – 18.43 | -27.89 | -84.70 – 28.91 | -24.82 | -82.69 – 33.05 | -22.19 | -80.19 – 35.81 | ||||||
| LAT CAT [Centro] | -31.00 *** | -40.16 – -21.83 | -26.66 *** | -35.84 – -17.48 | -9.65 | -25.97 – 6.66 | -8.56 | -24.85 – 7.72 | ||||||||
| LAT CAT [Sur] | -43.80 *** | -53.01 – -34.59 | -37.90 *** | -47.13 – -28.66 | -29.35 *** | -46.16 – -12.55 | -26.63 ** | -43.40 – -9.85 | ||||||||
| TIPO CARNADA [1] | -87.20 | -239.67 – 65.28 | ||||||||||||||
| TIPO CARNADA [2] | -50.33 | -202.64 – 101.98 | ||||||||||||||
| TIPO CARNADA [3] | 59.29 | -95.96 – 214.55 | ||||||||||||||
| zona [2] | -30.26 *** | -42.07 – -18.44 | ||||||||||||||
| zona [3] | -6.41 | -26.55 – 13.72 | ||||||||||||||
| prof med cat [501-1e+03] × zona2 |
6.93 | -18.62 – 32.47 | ||||||||||||||
| prof med cat [1001-1.5e+03] × zona2 |
-62.36 *** | -77.70 – -47.01 | ||||||||||||||
| prof med cat [1501-2e+03] × zona2 |
-10.01 | -27.22 – 7.20 | ||||||||||||||
| prof med cat [2001-3e+03] × zona2 |
-0.40 | -43.42 – 42.62 | ||||||||||||||
| prof med cat [501-1e+03] × zona3 |
24.17 | -7.32 – 55.66 | ||||||||||||||
| prof med cat [1001-1.5e+03] × zona3 |
-32.27 ** | -55.94 – -8.60 | ||||||||||||||
| prof med cat [1501-2e+03] × zona3 |
13.12 | -12.43 – 38.67 | ||||||||||||||
| prof med cat [2001-3e+03] × zona3 |
61.51 * | 6.00 – 117.02 | ||||||||||||||
| Observations | 3340 | 3340 | 3340 | 3340 | 3340 | 3293 | 3340 | 3340 | ||||||||
| R2 | 0.078 | 0.083 | 0.150 | 0.176 | 0.197 | 0.229 | 0.205 | 0.216 | ||||||||
| * p<0.05 ** p<0.01 *** p<0.001 | ||||||||||||||||
Table comparing performance model.
compare_performance(M01,
M02,
M03,
M04,
M05,
M06,
M07,
M08,
rank = TRUE,
verbose = FALSE)
# Comparison of Model Performance Indices
Name | Model | R2 | RMSE | Sigma | Performance-Score
----------------------------------------------------------
M06 | glm | 0.229 | 77.035 | 77.364 | 100.00%
M08 | glm | 0.216 | 77.282 | 77.667 | 94.44%
M07 | glm | 0.205 | 77.804 | 78.120 | 87.21%
M05 | glm | 0.197 | 78.206 | 78.500 | 81.47%
M04 | glm | 0.176 | 79.230 | 79.504 | 66.63%
M03 | glm | 0.150 | 80.476 | 80.670 | 48.76%
M02 | glm | 0.083 | 83.567 | 83.731 | 3.21%
M01 | glm | 0.078 | 83.796 | 83.921 | 0.00%Plot comparing performance model.
plot(compare_performance(M01,
M02,
M03,
M04,
M05,
M06,
M07,
M08,
rank = TRUE,
verbose = FALSE))model_performance(M06)
# Indices of model performance
R2 | RMSE | Sigma
-----------------------
0.229 | 77.035 | 77.364res_glm<-as.numeric(M06$residuals)
res_hist<-ggplot(as.data.frame(res_glm), aes(x=res_glm)) +
geom_histogram(position="identity", bins = 20)+
theme_bw()+
xlab("Residuos")+
ylab("Frecuencia")
#res_qqnorm<-qqnorm(res_glm)+
# qqline(res_glm)
res_density<-ggdensity(res_glm, res_glm = "", fill = "lightgray",
title = "") +
scale_x_continuous() +
stat_overlay_normal_density(color = "red", linetype = "dashed")+
xlab("Residuos")+
ylab("Densidad")
res_points<-ggplot(as.data.frame(res_glm),
aes(y=res_glm, x=(1:length(res_glm)))) +
geom_point()+
xlab("")+
ylab("Residuos")+
theme_few()
ggarrange(res_hist, res_density, res_points,
labels = c("A", "B", "C"),
ncol = 3)output<-M01$fitted.valuesRegarding this result, we can select best model and check some assumtion like:
- Colinealidad
- Data points influyentes
- Homoscedasticidad
- Normalidad de los residuos
- Independencia
COMPOSICIONES DE TALLAS
La información contenida en las estructiras de tallas son consideradas una de las piezas mas importantes en la ciencia pesquera (Canales, Punt, & Mardones, 2021; Hordyk, Ono, Sainsbury, Loneragan, & Prince, 2014; Rudd & Thorson, 2018). Sin embargo, esta debe ser validada y confirmada para no violar supuestos referidos a la repreentatividad del dato y su relación con la dinámica poblacional.
Ahora nos disponemos a explorar esta fuente de información.
Primero identificamos la estructura de la base;
glimpse(long)
Rows: 15,831
Columns: 24
$ ID <chr> "COD_BARCO=730881 AND TO_CHAR(FECHA_HORA_RECALADA,…
$ COD_BARCO <dbl> 730881, 730881, 730881, 730881, 730881, 730881, 73…
$ FECHA_HORA_RECALADA <dttm> 2021-05-27 15:00:00, 2021-05-27 15:00:00, 2021-05…
$ FECHA_HORA_ZARPE <dttm> 2021-05-10 11:00:00, 2021-05-10 11:00:00, 2021-05…
$ COD_PESQUERIA <dbl> 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50…
$ PUERTO_RECALADA <dbl> 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 3,…
$ NRO_FORMULARIO <dbl> 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9,…
$ REGION <dbl> 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1,…
$ FECHA_LANCE <dttm> 2021-05-27 15:00:00, 2021-05-27 15:00:00, 2021-05…
$ NUMERO_LANCE_EX <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
$ LATITUD <dbl> 243900, 243900, 243900, 243900, 243900, 243900, 24…
$ LONGITUD <dbl> 714400, 714400, 714400, 714400, 714400, 714400, 71…
$ NUMERO_EJEMPLARES <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ PESO_TOTAL_MUESTRA <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ NUMERO_CAJA <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
$ COD_ESPECIE <dbl> 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37…
$ ORIGEN_MUESTRA <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 3,…
$ FECHA_MUESTREO <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ N_TOTAL_INDIV <dbl> 85, 85, 85, 85, 85, 85, 85, 84, 84, 84, 84, 84, 23…
$ LONGITUD_DESCARTE <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ LONGITUD_MUESTRA <dbl> 115, 125, 129, 130, 143, 152, 157, 100, 108, 110, …
$ SEXO <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ N_INDIVIDUOS <dbl> 2, 1, 2, 2, 1, 1, 1, 2, 1, 2, 1, 3, 1, 1, 2, 1, 1,…
$ año <dbl> 2021, 2021, 2021, 2021, 2021, 2021, 2021, 2021, 20…
colSums(is.na(long))
ID COD_BARCO FECHA_HORA_RECALADA FECHA_HORA_ZARPE
2997 0 0 0
COD_PESQUERIA PUERTO_RECALADA NRO_FORMULARIO REGION
0 0 0 0
FECHA_LANCE NUMERO_LANCE_EX LATITUD LONGITUD
0 0 4425 4425
NUMERO_EJEMPLARES PESO_TOTAL_MUESTRA NUMERO_CAJA COD_ESPECIE
15733 15831 0 0
ORIGEN_MUESTRA FECHA_MUESTREO N_TOTAL_INDIV LONGITUD_DESCARTE
0 15831 0 15831
LONGITUD_MUESTRA SEXO N_INDIVIDUOS año
0 0 2 0Primero selecciono las columnas de interes y luego genero la expansión de LONGITUD_MUESTRA a N_INDIVIDUOS. Expand frecuency data related length, in this case N_INDIVIDUOS column have frecuency that we need expand to whole data frame.
Selecciono variables de interes
long1 <- long %>%
dplyr::select(6,8,11,12,21,22,23,24)
dim(long1)
[1] 15831 8long2 <- long %>%
drop_na(N_INDIVIDUOS) %>%
type.convert(as.is = TRUE) %>%
uncount(N_INDIVIDUOS)
dim(long2)
[1] 34299 23
table(long2$REGION)
1 3 4 5 7 8 10 12 14
494 3055 1697 12111 9452 4263 2848 1 378
legend.labels <- c('Indet.', 'Macho' , 'Hembra')
nbco <- ggplot(long2 %>%
filter(REGION!=12), aes(x=LONGITUD_MUESTRA, y = as.factor(año),
fill=as.factor(SEXO)))+
geom_density_ridges(stat = "density_ridges", bins = 20,
scale = 2, draw_baseline = FALSE,
alpha=0.5)+
facet_wrap(.~REGION, ncol=4) +
geom_vline(xintercept = 110, color = "red")+
scale_fill_viridis_d(option = "G",
name="SEXO",
labels=legend.labels)+
scale_y_discrete(breaks = seq(from = 2004, to = 2022, by = 2))+
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
theme_few()+
xlab("Longitud (cm.)")+
ylab("")
#scale_x_discrete((limits = rev(levels(talla2021$ANO_ARR))))+
nbcoPrepara los vectores para sumar a los .dat del modelo si los necesito.
long2$LONG_CAT <- as.numeric(as.character(cut(x = long2$LONGITUD_MUESTRA,
breaks = seq(10,220,2),
labels = seq(10,218,2),
right = FALSE)))
LONGT <- table(long2$año, long2$LONG_CAT)
# guardar en un .xlsBIOLOGICO
names(bio)
[1] "ID" "COD_BARCO"
[3] "FECHA_HORA_RECALADA" "FECHA_HORA_ZARPE"
[5] "PUERTO_RECALADA" "PUERTO_ZARPE"
[7] "REGION" "COD_PESQUERIA"
[9] "NRO_FORMULARIO" "FECHA_LANCE"
[11] "NUMERO_LANCE_EX" "LATITUD"
[13] "LONGITUD" "NUMERO_EJEMPLARES"
[15] "PESO_TOTAL_MUESTRA" "COD_ESPECIE"
[17] "N_ESPECIMEN" "ORIGEN_MUESTRA"
[19] "SEXO_ESPECIMEN" "NRO_MUESTRA"
[21] "LONGITUD_ESPECIMEN" "PESO_ESPECIMEN"
[23] "NUMERO_OTOLITOS" "PESO_EVISCERADO"
[25] "MADUREZ" "PESO_GONADAS"
[27] "CONTENIDO_ESTOMACAL" "REPLECION"
[29] "LONGITUD_DISCO" "ANCHO_DISCO"
[31] "PESO_ESTOMAGO" "PESO_PARED_ESTOMAGO"
[33] "NUMERO_FICHA" "CONT_ESTOMACAL_ESTADO_ESTOMAGO"
[35] "NUMERO_FICHA_GONADA" "año"
table(bio$COD_ESPECIE)
37
127409Filtro ciertos datos para estimar coeficientes
bio2 <- bio %>%
drop_na(PESO_ESPECIMEN,
LONGITUD_ESPECIMEN,
SEXO_ESPECIMEN) %>%
filter(PESO_ESPECIMEN>100,
año>2005)
Ahora por sexo
legend.label <- c('Macho' , 'Hembra')
my.plotsex<- ggplot(bio2 %>%
filter(SEXO_ESPECIMEN!="0",
SEXO_ESPECIMEN!="3",
SEXO_ESPECIMEN!="6"),
aes(y=PESO_ESPECIMEN, x=LONGITUD_ESPECIMEN,
colour=SEXO_ESPECIMEN)) +
geom_point(alpha=.1) +
scale_colour_manual(values = c("black", "red"),
name="Sexo",
labels=legend.label)+
facet_wrap(~ año, ncol = 4) +
geom_smooth(method = "loess",
se=TRUE,
formula = y ~ x,
size=0.3,
span=5)+
theme_few()+
theme(legend.position = "bottom")+
labs(y="Peso (gr.)",
x="Longitud (cm)")
my.plotsexAhora por región
legend.label <- c('Macho' , 'Hembra')
my.plotreg<- ggplot(bio2 %>%
drop_na(REGION) %>%
filter(SEXO_ESPECIMEN!="0",
SEXO_ESPECIMEN!="3",
SEXO_ESPECIMEN!="6",
REGION!=12),
aes(y=PESO_ESPECIMEN, x=LONGITUD_ESPECIMEN,
colour=SEXO_ESPECIMEN)) +
geom_point(alpha=.1) +
scale_colour_manual(values = c("black", "red"),
name="Sexo",
labels=legend.label)+
facet_wrap(~ REGION, ncol = 4) +
geom_smooth(method = "loess",
se=TRUE,
formula = y ~ x,
size=0.3,
span=5)+
theme_few()+
theme(legend.position = "bottom")+
labs(y="Peso (gr.)",
x="Longitud (cm)")
my.plotregExtraigo coeficientes globales
x<-log(bio2$LONGITUD_ESPECIMEN)
y<-log(bio2$PESO_ESPECIMEN)
reg.l<-lm(y~x)
summary(reg.l)
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-2.68627 -0.07255 -0.00070 0.07109 2.47970
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.005673 0.018393 -272.1 <2e-16 ***
x 3.082291 0.004136 745.3 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.1557 on 31985 degrees of freedom
Multiple R-squared: 0.9456, Adjusted R-squared: 0.9455
F-statistic: 5.555e+05 on 1 and 31985 DF, p-value: < 2.2e-16
anova(reg.l)
Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
x 1 13468.9 13469 555453 < 2.2e-16 ***
Residuals 31985 775.6 0
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
res <- residuals(reg.l)plot(res)
abline(h = 0, col = "red", lty=2)BITACORA FIPA 96- 14
bitfipa<- bit8696 %>%
mutate(CPUE = CAPTURA/DFP) %>%
drop_na(CPUE)
cpuefipa <- ggplot(bitfipa %>%
filter(ANO<1997) %>%
group_by(ANO) %>%
summarise(CPUEF=mean(CPUE)),
aes(ANO, CPUEF))+
geom_point(shape = 16, size = 3, alpha = 0.7)+
scale_x_continuous(breaks = seq(from = 1986, to = 2022, by = 1))+
geom_smooth(method = 'lm',
colour = 'blue',
size = 1.5)+
# stat_regline_equation(label.x=1993, label.y=500)+
# stat_cor(aes(label=..rr.label..), label.x=1993, label.y=450)+
#facet_wrap(~AREA)+
theme_bw()+
theme(axis.text.x = element_text(angle = 90, hjust = 2),
panel.grid = element_blank())+
labs(y="CPUE (kg/dfp)",
x="")
cpuefipaGuardo datos de la CPUE nominal del FIPA 96-14
cpuenomfipa <- bitfipa %>%
filter(ANO<1997) %>%
group_by(ANO) %>%
summarise(CPUEF=mean(CPUE))write_csv(cpuenomfipa, "CPUE_FIPA.csv")MAPAS
Lo primero es transformar los datos de lon_ctgr y lat_ctgr a formarto geom. A su vez, saco las NA de las coord
transformar los datos en un sf object
codmap <- st_as_sf(bitb2 %>%
drop_na(lon_ctgr) %>%
drop_na(lat_ctgr),
coords = c("lon_ctgr", "lat_ctgr"), crs = 4326)Ahora genero los mapas de Chile con un raster.
chile <- raster::getData("GADM", country = "CHL", level = 0)
chile1<-fortify(chile)
chilemap <- ggplot()+
geom_polygon(data=chile, aes(x=long, y=lat, group=group),
fill="lightblue",color="grey20", size=0.15)+
coord_sf(crs = st_crs(4326),
xlim = c(-80, -65),
ylim = c(-58, -15)) +
theme_bw()
chilemap# Aca veo los nombres
chile@data$NAME_1
NULLLuego genero los bordes sobre los cuales haré la grilla.
e <- extent(-80,-65,-58,-15)
rc <- crop(chile, e)
# la proyección adecuada en lat long
tcrs <- CRS("+init=epsg:4326")
transformed_points <- spTransform(rc, tcrs)
# para dejarlo en formato geom_sf
chile2 <- st_as_sf(transformed_points) Se estructura la grilla en el objeto raster de chile1
grid <- chile2 %>%
st_make_grid(cellsize = c(1,1)) %>% # para que quede cuadrada
st_cast("MULTIPOLYGON") %>%
st_sf() %>% # objeto en spatial feature
mutate(cellid = row_number())Pongo los datos codmap en la grilla. Aca solo elegí dfpy CPUE pero se ùeden resumir otros
joindat <- grid %>%
st_join(codmap) %>%
group_by(cellid, año) %>%
summarise(CPUEM = mean(CPUE),
DFPM =mean(dfp)) # por ejemplo, plotear la captura "CAPTURA_1"Mapa Esfuerzo
here we can viz diferent variables
## Plot final
mas1 <- ggplot() +
geom_sf(data=joindat %>%
filter(!is.na(DFPM)), aes(fill = DFPM),
color=NA) +
scale_fill_viridis_b(option="E",
direction=-1, name="DFP")+
geom_sf(data = grid, fill=NA, color=NA) +
geom_sf(data = chile2, color="grey", fill="white") +
coord_sf() +
geom_hline(yintercept = -47, color = "red")+
scale_alpha(guide="none")+
facet_wrap(~año, ncol=6)+
scale_x_continuous(breaks = seq(from = -80, to = -60, by = 10))+
xlab(expression(paste(Longitude^o,~'O'))) +
ylab(expression(paste(Latitude^o,~'S')))+
guides(colour = guide_legend()) +
theme_bw()+
theme(panel.background = element_rect(fill = 'white'),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
mas1Mapa CPUE
## Plot final
mas2 <- ggplot() +
geom_sf(data=joindat %>%
filter(!is.na(CPUEM)), aes(fill = CPUEM),
color=NA) +
scale_fill_viridis_b(option="G",
direction=-1, name="CPUE")+
geom_sf(data = grid, fill=NA, color=NA) +
geom_sf(data = chile2, color="grey", fill="white") +
coord_sf() +
geom_hline(yintercept = -47, color = "red")+
scale_alpha(guide="none")+
facet_wrap(~año, ncol=6)+
scale_x_continuous(breaks = seq(from = -80, to = -60, by = 10))+
xlab(expression(paste(Longitude^o,~'O'))) +
ylab(expression(paste(Latitude^o,~'S')))+
guides(colour = guide_legend()) +
theme_bw()+
theme(panel.background = element_rect(fill = 'white'),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
mas2Falta componer mapas de tallas medias y alguna otra variable de interés.
En función de los datos analizados, el modelo propuesto es un enfoque de producción global estado espacio descrito por Pedersen & Berg (2017) y Acom (2015).